%this functions calculates the velocity ensemble average of several
%experiments. The inputs are:
%direction: main flow direction w.r.t the CCD
%pathFiles: vector with the entire path of the vector fields
%pathJunk: vector with the entire path of the wrong vector fields
%windowProcess: sub-region of the field to be processed

%created by Ricardo Mejia-Alvarez. Urbana, IL. 05/06/08.

%modified:  07/10/08
%           02/14/09
%           02/20/09
%           08/01/09
%           12/20/09
%           01/02/2010
%           01/25/2010
%           04/06/2010
%           04/18/2010
%           05/30/2010


function [Xc,Yc,Zc,Uxavrg,Uyavrg,Uzavrg,CHCtotal,I,J,nc,FileNames] = ensembleParallel(direction,pathFiles,pathJunk,windowProcess,FAC,invertCoord)
N = 0;      %counter for the number of files used in the ensemble average

thr1 = 10;   %threshold to filter wrong fields based on standarized
            %residuals of each velocity component within instantaneous
            %realizations

thr2 = 10;   %threshold to filter wrong fields based on standarized
            %residuals of squares of each velocity component within
            %instantaneous realizations

thr3 = 10;   %threshold to filter wrong fields based on standarized
            %residuals of the differences between the maximum and median
            %values of the squared values of velocity components. This
            %filtering is carried out taking the whole ensamble into
            %account
            

path = char(pathFiles);


UU(1:length(path(:,1)) , 3) = 0;
UUmax(1:length(path(:,1)) , 3) = 0;
FileNames(1:length(path(:,1))) = {0};


    for j = 1 : length(path(:,1))
        
        
        if j > 1
            
            if nc == 5
                [nc,I,J,~,~,~,~,Ux,Uy,CHC] = matrix( path(j,:) );
                [I,J,~,~,Ux,Uy,CHC] = WindowFile(windowProcess,I,J,Xc,Yc,Ux,Uy,CHC);
                if FAC ~= 1
                    Ux = Ux * FAC;
                    Uy = Uy * FAC;
                end
                Uz = zeros(J,I);
                
                Ux = Ux * invertCoord(1);
                Uy = Uy * invertCoord(2);
                
            elseif nc == 4
                [~,~,~,~,~,~,~,Ux,Uy] = matrix( path(j,:) );
                [~,~,~,~,Ux,Uy,CHC] = WindowFile(windowProcess,I,J,Xc,Yc,Ux,Uy,CHC);
                if FAC ~= 1
                    Ux = Ux * FAC;
                    Uy = Uy * FAC;
                end
                Uz = zeros(J,I);
                
                Ux = Ux * invertCoord(1);
                Uy = Uy * invertCoord(2);
                
            else
                [~,~,~,~,~,~,~,~,Ux,Uy,Uz,CHC] = matrix( path(j,:) );
                [~,~,~,~,~,Ux,Uy,Uz,CHC] = WindowFile(windowProcess,I,J,Xc,Yc,Zc,Ux,Uy,Uz,CHC);
                if FAC ~= 1
                    Ux = Ux * FAC;
                    Uy = Uy * FAC;
                    Uz = Uz * FAC;
                end
                
                Ux = Ux * invertCoord(1);
                Uy = Uy * invertCoord(2);
                Uz = Uz * invertCoord(3);
            end
            
            [CHC,rowReplace,colReplace] = depurationCHC(CHC);   %identifies places with bad vectors as generated by Insight 3G flag
            
            for k = 1 : length(rowReplace)
                Ux(rowReplace(k) , colReplace(k)) = 0;
                Uy(rowReplace(k) , colReplace(k)) = 0;
                Uz(rowReplace(k) , colReplace(k)) = 0;
            end
            
            %any place with a wrong vector has now a 0 and is not going to
            %participate in the ensemble averaging.
            Ux = Ux .* CHC;
            Uy = Uy .* CHC;
            Uz = Uz .* CHC;
            
            if direction ~= 3
                if direction == 1
                    U1 = Ux;
                    U2 = Uy;
                    U3 = Uz;
                    U4 = Ux.^2;
                    U5 = Uy.^2;
                    U6 = Uz.^2;
                else
                    U1 = Ux';
                    U2 = Uy';
                    U3 = Uz';
                    U4 = (Ux').^2;
                    U5 = (Uy').^2;
                    U6 = (Uz').^2;
                end
                
                [wrong] = depurationResiduals(thr1,path(j,:),pathJunk,U1,U2,U3);
                if wrong == 0
                    [wrong] = depurationResiduals(thr2,path(j,:),pathJunk,U4,U5,U6);
                end
                
            else
                wrong = 0;
            end
                
            
            if wrong == 0
                Uxavrg = Uxavrg + Ux;
                Uyavrg = Uyavrg + Uy;
                Uzavrg = Uzavrg + Uz;
                
                N = N + 1;
                CHCtotal = CHCtotal + CHC;
                FileNames(N) = {path(j,:)};
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %calculating the medians and maxima of the square values of
                %instantaneous velocity components fields. This will be
                %used to exclude fields with strong bias
                [UU(N,:) , UUmax(N,:)] = medianDep(Ux,Uy,Uz);
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
            else
            end
              
 
        elseif j == 1
            
            [nc] = matrix( path(j,:) );
            
            if nc == 5
                [nc,I,J,~,~,Xc,Yc,Uxavrg,Uyavrg,CHC] = matrix( path(j,:) );
                [I,J,Xc,Yc,Uxavrg,Uyavrg,CHC] = WindowFile(windowProcess,I,J,Xc,Yc,Uxavrg,Uyavrg,CHC);
                Zc = zeros(J,I);
                if FAC ~= 1
                    Uxavrg = Uxavrg * FAC;
                    Uyavrg = Uyavrg * FAC;
                end
                Uzavrg = zeros(J,I);
                
                Xc = Xc * invertCoord(1);
                Yc = Yc * invertCoord(2);
                Uxavrg = Uxavrg * invertCoord(1);
                Uyavrg = Uyavrg * invertCoord(2);

            elseif nc == 4
                [nc,I,J,~,~,Xc,Yc,Uxavrg,Uyavrg] = matrix( path(j,:) );
                CHC = Xc ./ Xc;
                [I,J,Xc,Yc,Uxavrg,Uyavrg,CHC] = WindowFile(windowProcess,I,J,Xc,Yc,Uxavrg,Uyavrg,CHC);
                Zc = zeros(J,I);
                if FAC ~= 1
                    Uxavrg = Uxavrg * FAC;
                    Uyavrg = Uyavrg * FAC;
                end
                Uzavrg = zeros(J,I);
                
                Xc = Xc * invertCoord(1);
                Yc = Yc * invertCoord(2);
                Uxavrg = Uxavrg * invertCoord(1);
                Uyavrg = Uyavrg * invertCoord(2);
                
            else
                [nc,I,J,~,~,Xc,Yc,Zc,Uxavrg,Uyavrg,Uzavrg,CHC] = matrix( path(j,:) );
                [I,J,Xc,Yc,Zc,Uxavrg,Uyavrg,Uzavrg,CHC] = WindowFile(windowProcess,I,J,Xc,Yc,Zc,Uxavrg,Uyavrg,Uzavrg,CHC);
                if FAC ~= 1
                    Uxavrg = Uxavrg * FAC;
                    Uyavrg = Uyavrg * FAC;
                    Uzavrg = Uzavrg * FAC;
                end
                Xc = Xc * invertCoord(1);
                Yc = Yc * invertCoord(2);
                Zc = Zc * invertCoord(3);
                Uxavrg = Uxavrg * invertCoord(1);
                Uyavrg = Uyavrg * invertCoord(2);
                Uzavrg = Uzavrg * invertCoord(3);
            end
                       
            [CHC,rowReplace,colReplace] = depurationCHC(CHC);   %identifies places with bad vectors as generated by Insight 3G flag
            
            for k = 1 : length(rowReplace)
                Uxavrg(rowReplace(k) , colReplace(k)) = 0;
                Uyavrg(rowReplace(k) , colReplace(k)) = 0;
                Uzavrg(rowReplace(k) , colReplace(k)) = 0;
            end
            
            %any place with a wrong vector has now a 0 and is not going to
            %participate in the ensemble averaging.
            Uxavrg = Uxavrg .* CHC;
            Uyavrg = Uyavrg .* CHC;
            Uzavrg = Uzavrg .* CHC;

            if direction ~= 3
                if direction == 1
                    U1 = Uxavrg;
                    U2 = Uyavrg;
                    U3 = Uzavrg;
                    U4 = Uxavrg.^2;
                    U5 = Uyavrg.^2;
                    U6 = Uzavrg.^2;
                else
                    U1 = Uxavrg';
                    U2 = Uyavrg';
                    U3 = Uzavrg';
                    U4 = (Uxavrg').^2;
                    U5 = (Uyavrg').^2;
                    U6 = (Uzavrg').^2;
                end
                
                [wrong] = depurationResiduals(thr1,path(j,:),pathJunk,U1,U2,U3);
                if wrong == 0
                    [wrong] = depurationResiduals(thr2,path(j,:),pathJunk,U4,U5,U6);
                end
                                
            else
                wrong = 0;
            end
            

            if wrong == 0
                CHCtotal = CHC;
                N = N + 1;
                FileNames(N) = {path(j,:)};
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %calculating the medians and maxima of the square values of
                %instantaneous velocity components fields. This will be
                %used to exclude fields with strong bias
                [UU(N,:) , UUmax(N,:)] = medianDep(Uxavrg,Uyavrg,Uzavrg);
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
            else
                Uxavrg = zeros(J,I);
                Uyavrg = zeros(J,I);
                Uzavrg = zeros(J,I);
                CHCtotal = zeros(J,I);
            end
 
        end
        
    end
    
    FileNames = FileNames(1:N);
    FileNames = char(FileNames);
    
    UU = UU(1:N , :);
    UUmax = UUmax(1:N , :);
    
    %filtering the data based on the differences between the maximum
    %instantaneous values of UUi and their instantaneous medians
    
    ResUU = UUmax - UU;
    iqrResUU = iqr(ResUU);
    iqrResUU = repmat(iqrResUU , [length(ResUU(:,1)) , 1]); 
    NormRes = ResUU ./ iqrResUU;
    
    %find the indices of the files that do not satisfy the criteria to be
    %accepted.
    [row , col] = find( abs( NormRes ) > thr3); %#ok<NASGU>
    clear col
    test = unique(row);
    
    
    if isempty(test) == 1
    else
        
        for k = 1 : length(test)
            fid = test(k);
            badFile = FileNames(fid,:);
            display(badFile)
            
            if nc == 5
                [nc,I,J,~,~,Xc,Yc,Ux,Uy,CHC] = matrix( badFile );
                [I,J,Xc,Yc,Ux,Uy,CHC] = WindowFile(windowProcess,I,J,Xc,Yc,Ux,Uy,CHC);
                Zc = zeros(J,I);
                Uz = zeros(J,I);
            else
                [nc,I,J,~,~,Xc,Yc,Zc,Ux,Uy,Uz,CHC] = matrix( badFile );
                [I,J,Xc,Yc,Zc,Ux,Uy,Uz,CHC] = WindowFile(windowProcess,I,J,Xc,Yc,Zc,Ux,Uy,Uz,CHC);
            end
            
            [CHC,rowReplace,colReplace] = depurationCHC(CHC);   %identifies places with bad vectors as generated by Insight 3G flag
            
            for m = 1 : length(rowReplace)
                Ux(rowReplace(m) , colReplace(m)) = 0;
                Uy(rowReplace(m) , colReplace(m)) = 0;
                Uz(rowReplace(m) , colReplace(m)) = 0;
            end
            
            Ux = Ux .* CHC;
            Uy = Uy .* CHC;
            Uz = Uz .* CHC;
            
      
            Uxavrg = Uxavrg - Ux;
            Uyavrg = Uyavrg - Uy;
            Uzavrg = Uzavrg - Uz;
            CHCtotal = CHCtotal - CHC;
            
            movefile(badFile,pathJunk);
        end

    end
    
   
    Uxavrg = Uxavrg ./ CHCtotal;
    Uyavrg = Uyavrg ./ CHCtotal;
    Uzavrg = Uzavrg ./ CHCtotal;
    
    [~,rowReplace,colReplace] = depurationCHC(CHCtotal);
    
    for m = 1 : length(rowReplace)
                Uxavrg(rowReplace(m) , colReplace(m)) = 0;
                Uyavrg(rowReplace(m) , colReplace(m)) = 0;
                Uzavrg(rowReplace(m) , colReplace(m)) = 0;
    end
    
    
    FileNames = char(FileNames);
    firstFileName = filename(FileNames(1,:));
    FileNames = FileNames(: , length(FileNames(1,:)) - length(firstFileName) + 1 : length(FileNames(1,:)) );
    
end